1. Field of the Invention
The present invention relates to an interpolation method that is used in computational fluid dynamics for realizing the motion of a fluid on a computer, an apparatus for carrying out the method, and a program for implementing the method.
2. Description of the Related Art
In recent years, computational fluid dynamics that realizes the motion of a fluid on a computer has made remarkable progress. However, it is difficult to describe the motion of an object in a fluid, or the motion of two fluids or a multi-layer fluid, and therefore, various approaches have been made to this difficult theme.
To compute a fluid flow around a stationary object, a grid can be generated according to the shape of the object, and it is a conventionally known method of the computation to utilize an unstructured grid or the like based on a finite element method or the like. When the object is moved, a method of changing the grid according to the motion of the object can be employed. However, if the change of the grid involves large deformation, it is generally difficult to perform the computation.
On the other hand, a difference equation that approximately realizes a fluid equation describing a fluid flow on a structured grid, by using duality between a grid (primary grid) and a dual grid known as a staggered grid, is highly reliable due to its stability. This leads to an idea of describing a volume of a substance in each cell as a function on a grid, and moving the function. One example of the method based on this idea utilizes the Volume-of-Fluid method (hereinafter referred to as “the VOF method”) for computing a flow of an incompressible fluid. In this case, it is a key technique to move a function describing the volume of the substance of the fluid.
Further, in the case of a compressible fluid, the computation of a flow thereof involves key problems of defining motions of various objects, and in particular, propagation of a shock wave and the like are very closely related with the above-mentioned theme of describing the motion of a substance or an interface.
An initial-value problem concerning an equation (referred to as “an advection equation”) describing such a motion is called the Riemann problem after a paper by Riemann on a wave equation concerning a wave in the air “Ueber die Fortpflanzung ebener Luftwellen von endlicher Schwingungsweite (Albhandlungen der Koniglichen Gesellschaft der Wissenschaften zu Gottingen, 8 (1860), 43-65”.
As numerical solutions to the Riemann problem, there have been proposed various concepts and schemes, independent of each other or dependent on each other, including a high-order accuracy upwind scheme called the Total Variation Diminishing scheme (hereinafter referred to as “the TVD scheme”), the Essentially Non-Oscillatory scheme (hereinafter referred to as “the ENO scheme”) which introduces a limiter function, the Monotone Upstream Centered Scheme for Conservation Laws (hereinafter referred to as “the MUSCL scheme”, and the Positivity-Preserving scheme (hereinafter referred to as “the PP scheme”.
One of equations dealt with in the present invention is a one-dimensional partial differential equation defined by:∂F(x,t)/∂t=∂C(x,t)F(x,t)/∂x  (1)
To perform computation based on the equation, we will now consider the problem of converting this equation into a difference equation. In the following, arguments X, T of a function F(X,T) represent a discrete grid point in space and a discrete time point, respectively.
To solve the equation (1), there have been proposed various computing methods in the computational fluid dynamics as described above. The TVD scheme is to perform the computation under a TVD condition that a TVD value TV(F[•,n]) defined by:TV(F[•,n])=Σj|F[j+1,n]−F[j,n]|  (2)does not increase with time evolution, the TVD condition being defined by the following inequality:TV(F[•,n+1])≦TV(F[•,n])  (3)
The MUSCL scheme, in a broader sense, is a computational method composed of data reconstruction and an approximate Riemann solution.
The ENO scheme is an upwind method with a limiter function in which the MUSCL scheme and the TVD scheme are combined to eliminate vibration occurring when there are discontinuities in the above-described function.
Further, for data reconstruction, there are conventionally employed a method of interpolating the above-described function with a piecewise constant function (which gives first-order accuracy), a method of interpolating the above-described function with a piecewise linear function (which gives second-order accuracy), and a method of interpolating the above-described function by a piecewise parabolic function (which gives third-order accuracy).
The present invention relates to improvement in the second-order accuracy of interpolation.
Now, a description will be given of a conventional second-order ENO scheme.
First, a description will be given of a method of numerically solving the above equation (1) by finding a solution thereto by finite difference approximation, assuming that there is a piecewise linear function:H[I,J](x)=F(I,J)+G(I,J)x  (4)
If the function C(I+½, J) is positive, there is given the following definition:Q(I+½,J):=C(I+½,J)H[I,J](Δx/2−C(I,J)Δt/2)  (5)
If the function is not positive, there is given the following definition:
                              Q          ⁢                                          ⁢                      (                                          I                +                                  1                  /                  2                                            ,              J                        )                          :=                  C          ⁢                                          ⁢                      (                                          I                +                                  1                  /                  2                                            ,              J                        )                    ⁢                                          ⁢                      H            ⁡                          [                                                I                  +                  1                                ,                J                            ]                                ⁢                                          ⁢                      (                                                            -                  Δ                                ⁢                                                                  ⁢                                  x                  /                  2                                            -                                                C                  ⁡                                      (                                                                  I                        +                        1                                            ,                      J                                        )                                                  ⁢                                                                  ⁢                Δ                ⁢                                                                  ⁢                                  t                  /                  2                                                      )                                              (        6        )            
The above equation (1) can be approximately solved by obtaining, from the above definitions, the following equation:
                              F          ⁢                                          ⁢                      (                          I              ,                              J                +                1                                      )                          =                              F            ⁢                                                  ⁢                          (                              I                ,                J                            )                                +                      Δ            ⁢                                                  ⁢            t            ⁢                                                  ⁢                                          (                                                      Q                    ⁢                                                                                  ⁢                                          (                                                                        I                          +                                                      1                            /                            2                                                                          ,                        J                                            )                                                        -                                      Q                    ⁢                                                                                  ⁢                                          (                                                                        I                          -                                                      1                            /                            2                                                                          ,                        J                                            )                                                                      )                            /              Δ                        ⁢                                                  ⁢            x                                              (        7        )            
On the other hand, a slope function G(I,J) representative of the slope of the piecewise linear function of the equation (4) is defined by the following equation:G(I,J)=ΔF(I,J)/Δx  (8)
Now, there are various methods of defining ΔF(I,J), as described hereinafter, and typical three of them will be described here.
First, there are introduced the following definitions:Sp(I,J):=F(I+1,J)−F(I,J)Sm(I,J):=F(I,J)−F(I−1,J)  (9)
A TVD limiter function in one dimension called the minmod function gives ΔF(I,J) defined by:ΔF(I,J)=(sgn(Sp(I,J))+sgn(Sm(I,J)))Min(|Sp(I,J)|,|Sm(I,J)|)/2  (10)
A TVD limiter function in one dimension called the averaging limiter function gives ΔF(I,J) defined by:ΔF(I,J)=(sgn(Sp(I,J))+sgn(Sm(I,J)))Min(|Sp(I,J)+Sm(I,J)|/2,2|Sp(I,J)|,2|Sm(I,J)|)/2  (11)
A TVD limiter function in one dimension called the superbee limiter function gives ΔF(I,J) defined by:ΔF(I,J)=(sgn(Sp(I,J))+sgn(Sm(I,J)))Min(Max(|Sp(I,J)|,|Sm(I,J)|)),2|Sp(I,J)|,2|Sm(I,J)|)/2  (12)
In the above equations (10) to (12), sgn( ) represents a sign function. That is, when the argument thereof is positive i.e. of a plus sign, the function gives a value of 1, while the same is negative i.e. of a minus sign, the function gives a value of −1.
These methods are described e.g. in “A. Suresh, “Positivity-Preserving Schemes in Multidimensions” SIAM J. Sci. Comput. 22, 1184-1198(2000)”.
Using these interpolation methods, a function describing a square shape was subjected to advection calculation using the above equation (7). FIGS. 7 to 9 show results of the calculation. More specifically, in FIGS. 7 to 9, shapes are depicted which are given by the function after 100 steps under periodic boundary conditions, using the minmod limiter function, the averaging limiter function, and the superbee limiter function, together with the initial shape for comparison. As is clear from FIGS. 7 to 9, the square shape as the initial shape is drastically lost into dull shapes.
Next, a description will be given of an algorithm for a higher dimension obtained by extending the algorithm for one dimension described heretofore. A method employed in the algorithm for a higher dimension described here is disclosed in the above-mentioned literature.
Suppose that F(I,J,T) represents a function of a point on a two-dimensional grid at a given time T, wherein I and J are integers and cooperatively represent a grid point on the two-dimensional grid.
For the function, there are defined the following values:Vmin=Min(F(I−1,J+1,T)−F(I,J,T),F(I,J+1,T)−F(I,J,T),F(I+1,J+1,T)−F(I,J,T),F(I−1,J,T)−F(I,J,T),−ε, F(I+1,J,T)−F(I,J,T),F(I−1,J−1,T)−F(I,J,T),F(I,J−1,T)−F(I,J,T),F(I+1,J−1,T)−F(I,J,T))Vmax=Max(F(I−1,J+1,T)−F(I,J,T),F(I,J+1,T)−F(I,J,T),F(I+1,J+1,T)−F(I,J,T),F(I−1,J,T)−F(I,J,T),ε,F(I+1,J,T)−F(I,J,T),F(I−1,J−1,T)−F(I,J,T),F(I,J−1,T)−F(I,J,T),F(I+1,J−1,T)−F(I,J,T))  (13)
Further, there are defined:Sx(I,J,T)=(F(I+1,J,T)−F(I−1,J,T))/2Sy(I,J,T)=(F(I,J+1,T)−F(I,J−1,T))/2  (14)V=2(Min(|Vmin|,|Vmax|))/(|Sx(I,J,T)|+|Sx(I,J,T)|+ε)  (15)
At this time, after defining:Gx(I,J,T)=Min(1,V)Sx(I,J,T)/ΔxGy(I,J,T)=Min(1,V)Sy(I,J,T)/Δy  (16)an interpolation function H(I,J,T)(x,y) for a grid point (I,J) is defined as:H(I,J,T)(x,y)=F(I,J,T)+Gx(I,J,T)x+Gy(I,J,T)y  (17)
By using the method, with respect to initial values of a function shown in FIG. 10A, advection calculation of the function can be carried, a result of which is shown in FIG. 10B.
However, the conventional interpolation method suffers from the following problem:
In the Volume-of-Fluid (VOF) method applied to a two-phase incompressible fluid, how to cause time evolution of a shape-describing function descriptive of a shape as a result of the advection calculation while preserving the shape described by the function is a matter of great interest, but the conventional interpolation method tends to cause progressive loss of the initial shape into a dull shape, with execution of computational steps.
In view of this problem, a computational algorithm enabling the shape to be preserved more effectively is indispensable, and to this end, for an algorithm for one dimension, it is necessary to replace the interpolation function of the equation (4) by an improved one, and for an algorithm for a higher dimension, it is necessary to replace the interpolation function of the equation (17) by an improved one.